##### WTP analysis smartphones ######


rm(list = ls())


# Load required packages
library(cregg)
library(ggplot2)
library(dfidx)
library(mlogit)
library(dplyr)
library(readr)
library(stargazer)
library(performance)

# Load conjoint data
dfconj <- read_csv("/Users/cbrugge/Desktop/CE Paper/wtp_data_ce.csv") %>% 
  as.data.frame() %>% 
  select(-"...1")


# Prep data for dfidx  ----------------------------------------------------

dfconj_prod1 <- subset(dfconj, product == 1)

colnames(dfconj_prod1) <- c("ID", "product" , "roundid", "choice", "cost.A", "life.A", "env.A", 
                           "eff.A", "recy.A", "rep.A", "product.type")

dfconj_prod1 <- subset(dfconj_prod1, select = -c(product))

dfconj_prod2 <- subset(dfconj, product == 2)

colnames(dfconj_prod2) <- c("ID", "product" , "roundid", "choice", "cost.B", "life.B", "env.B", 
                           "eff.B", "recy.B", "rep.B", "product.type")


dfconj_prod2 <- subset(dfconj_prod2, select = -c(product, choice, roundid, ID, product.type))

data <- cbind(dfconj_prod1, dfconj_prod2)

col_order <- c("ID", "roundid", "choice", "cost.A", "life.A", "env.A", "eff.A", "recy.A", "rep.A", 
               "cost.B", "life.B", "env.B", "eff.B", "recy.B", "rep.B", "product.type")

data <- data[, col_order]

# convert relevant variables to characters
data <- data %>%
  mutate(ID = as.character(ID),
         roundid = as.character(roundid),
         choice = as.character(choice))

#add character string to ID 
data$ID <- paste("R_", data$ID, sep = "")

data <- data %>% mutate(choice = ifelse(choice == "1", "A", "B"))

# create individual datasets by product -----------------------------------
data_smartphone <- subset(data, product.type == "Smartphone")

# Create a new unique identifier by combining "round_id" and "ID"
data_smartphone$unique_id <- as.character(interaction(data_smartphone$roundid, data_smartphone$ID, sep = "."))


# apply dfidx --------------------------------------------------------------

mlogdata <- dfidx(data_smartphone, idx = "unique_id", choice = "choice", varying = 4:15, sep = ".") # dfidx takes as first argument a data.frame and returns a dfidx object, which is a data.frame in "long" format with a special data frame column which contains the indexes


# Convert them to the new values 250, 850, 1900, and 3100 to allow for non-linear analysis
mlogdata <- mlogdata %>%
  mutate(cost = case_when(
    cost == 1 ~ 180,
    cost == 2 ~ 550,
    cost == 3 ~ 920 ,
    cost == 4 ~ 1300,
    TRUE ~ cost  # keep unchanged if not matching any condition
  ))
# calculate mlogit and wtp ------------------------------------------------

# general model

m1 <- mlogit(choice ~ cost + life + env + eff + recy + rep | -1, 
             mlogdata)

summary(m1)

# Print & export regression table
stargazer(m1, type = "text", covariate.labels = c("Price","Function duration","Environmental impact during production",
          "Energy efficiency","Recyclability","Repairability"), dep.var.labels =
            "Product choice")


stargazer(m1, type = "html", covariate.labels = c("Price","Function duration","Environmental impact during production",
                                                  "Energy efficiency","Recyclability","Repairability"), dep.var.labels =
            "Product choice (smartphone)", out = "/Users/cbrugge/Desktop/CE Paper/mixed_logit_smartphones.html")

# odds ratios
odds <- exp(coef(m1, matrix = TRUE))
precentage_increase_by_level <- 1-odds
precentage_increase_by_level

# calculate wtp

wtp <- coef(m1)[2:length(coef(m1))] / coef(m1)["cost"]
#After running this code, the wtp variable will contain the calculated WTP values
#for each attribute level (excluding the intercept) relative to a one-point increase
#in the price attribute.

# Transform the WTP values to percentage change
wtp_percentage_change <- (exp(wtp) - 1) * 100
print(wtp_percentage_change)


